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Abstract.  Recently  some  methods  have  been  presented  to  extract 
ordinary  differential  equations  (ODE)  directly  from  an  experimental 
time  series.  Here,  we  introduce  a  new  method  to  find  an  ODE  which 
models  both  the  short  time  and  the  long  time  dynamics.  The  exper¬ 
imental  data  are  represented  in  a  state  space  and  the  corresponding 
flow  vectors  are  approximated  by  polynomials  of  the  state  vector  com¬ 
ponents.  We  apply  these  methods  both  to  simulated  data  and  experi¬ 
mental  data  from  human  limb  movements,  which  like  many  other  bio¬ 
logical  systems  can  exhibit  limit  cycle  dynamics.  In  systems  with  only 
one  oscillator  there  is  excellent  agreement  between  the  limit  cycling 
displayed  by  the  experimental  system  and  the  reconstructed  model, 
even  if  the  data  are  very  noisy.  Furthermore  we  study  systems  of 
two  coupled  limit  cycle  oscillators.  There,  a  reconstruction  was  only 
successful  for  data  with  a  sufficently  long  transient  trajectory  and  rel¬ 
atively  low  noise  level. 

1.  Introduction 

The  modeling  and  analysis  of  dynamical  systems  is  a  field  of  increasing  inter¬ 
est,  in  part  because  of  applications  in  forecasting  and  control  [?].  For  a  long 
time  linear  models  have  dominated  descriptions  of  dynamical  systems  and 
control  theoretic  approaches  [?].  Very  complex  and  randomlike  behavior  was 
viewed  from  a  statistical  perspective  in  which  very  many  degrees  of  freedom 
were  involved.  Only  quite  recently  have  nonlinear  models  emerged,  capable 
of  mirroring  chaotic  dynamics  and  other  phenomena  such  as  self  excited  oscil¬ 
lation  [?].  Such  systems  can  exhibit  extremely  complex  dynamical  behavior, 
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even  though  the  underlying  dynamics  may  be  low  dimensional.  On  the  other 
hand  very  high  dimensional  systems  such  as  fluids  or  lasers  can  show  simple 
and  low  dimensional  dynamics,  which  may  be  described  using  low  dimen¬ 
sional  models  [4,5].  In  such  systems  most  degrees  of  freedom  are  slaved  and 
only  a  small  number  of  order  parameters  are  necessary  for  a  description  [5]. 

Along  with  new  theoretical  concepts  have  come  a  variety  of  different  tech¬ 
niques  to  characterize  the  dynamics  of  such  complex  systems  from  any  time 
series  representative  of  the  system’s  dynamical  behavior.  For  example,  gen¬ 
eralized  dimensions  or  Liapunov  exponents  have  been  used  and  calculations 
performed  in  complex  systems  all  the  way  from  the  laser  to  the  human  EEG 
[6]. 

In  this  article,  we  present  a  new  method  for  constructing  a  model  equation 
from  experimental  time  series  that  reproduces  the  dynamical  behavior  of  the 
given  system.  The  time  series  has  to  be  represented  in  a  low  dimensional 
state  space,  e.g.  position  and  velocity  for  a  one  dimensional  oscillator.  A 
trajectory  of  the  system  is  given  by  the  evolution  of  a  point  in  the  state 
space.  A  flow  vector  field,  which  is  a  function  of  the  state  vector  components, 
governs  this  dynamical  behavior.  In  order  to  model  the  given  time  series  with 
a  low  dimensional  ordinary  differential  equation  (ODE)  or  flow  vector,  the 
dynamics  has  to  be  deterministic,  i.e.  no  crossovers  of  the  trajectory  in  the 
state  space  exist  that  are  not  due  to  noise  [7].  Furthermore,  the  mapping  of 
the  given  time  series  into  state  space  has  to  be  bijective.  That  is,  a  discrete 
point  to  point  mapping  must  exist  between  the  trajectory  given  as  a  time 
series,  and  the  trajectory  given  in  state  space,  and  vice  versa.  In  one  of  our 
experimental  examples  from  human  limb  movement,  pauses  occur  which  do 
not  follow  this  condition. 

The  reconstructed  model  equation,  e.g.  an  ODE,  is  a  “global  description” 

[8] ,  meaning  the  equation  is  valid  on  the  whole  subspace  of  the  state  space 
where  experimental  data  are  available.  The  model  equation  is  represented  in 
polynomials  of  the  components  of  the  state  vector  and  the  coefficients  of  these 
polynomials  are  fitted  to  the  experimental  time  series.  Differential  equations 
of  this  t)’pe  are  always  deterministic  since  they  obey  the  Lipschitz  condition 

[9] .  An  advantage  of  this  closed  model  equation  is  that  it  can  be  compared  to 
theoretical  models  and  used  to  analyze  and  interpret  the  system’s  behavior. 
Furthermore,  this  description  may  be  implemented  to  control  a  given  system 

[10] . 

A  disadvantage  of  a  global  description  is  that  for  higher  dimensional  sys¬ 
tems  and  for  higher  order  fits,  the  number  of  coefficients  increases  rapidly, 
rendering  the  approach  cumbersome  and  unfavorable.  In  this  regard  we  want 
to  mention  another  approach,  in  which  the  dynamics  is  represented  by  a  “lo¬ 
cal  description”  [2,8,11]  Here  the  trajectory  is  forecast  by  estimating  the 
future  trajectory  by  the  trajectories  of  past  states  that  come  close  to  the 
present  state.  Thus,  the  model  equations  are  only  valid  in  a  small  area  of 
state  space,  i.e.,  no  closed  description  is  possible.  On  the  other  hand,  it 
may  be  possible  to  forecast  higher  dimensional  systems  better  than  using  the 
global  description. 
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The  goal  of  the  present  work  is  to  provide  a  global  description  for  a 
given  experimental  time  series,  by  reconstructing  the  underlying  dynamics 
eis  well  as  possible.  Of  course  the  resulting  model,  in  general,  does  not 
provide  an  explanation  for  the  system’s  dynamics.  Nevertheless,  at  least 
in  simple  cases,  it  may  still  be  possible  to  connect  the  dynamics  and  the 
model  equation  to  a  deeper  level  of  description,  one  which  would  allow  for 
an  understanding  of  the  reconstructed  equation.  As  an  aside,  such  insight 
can  be  enhanced  by  empirical  study.  For  example  in  [12]  we  show  how  an 
observable,  the  relaxation  time  (i.e.  the  time  taken  to  return  to  the  limit  cycle 
after  a  perturbation),  is  related  to  the  strength  of  the  system’s  nonlinearity. 
In  our  approach  only  very  obvious  state  vector  components  are  used,  namely 
the  input  data  and  their  derivatives.  Other  state  vectors,  such  as  those  that 
can  be  created  using  delays  [13],  are  less  helpful  in  analyzing  and  interpreting 
the  reconstructed  equations. 

In  earlier  work  the  model  equation  was  obtained  by  a  direct  fit  of  the 
flow  vector  field,  using  an  approach  developed  by  Cremers  and  Hiibler  [7]. 
A  similar  method  was  investigated  in  a  slightly  different  manner  by  Crutch¬ 
field  and  McNamara  [8].  The  Cremers  and  Hiibler  approach,  from  now  on 
referred  to  as  the  “flow  method”,  produces  very  good  results  if  data  from 
a  large  enough  subspace  of  the  statespace  are  available.  However,  in  many 
experimental  applications  the  system  is  not  chaotic;  rather  it  is  attracted 
very  rapidly  to  a  limit  cycle.  In  such  cases,  only  data  from  the  limit  cycle 
can  be  obtained,  in  the  complete  or  near  absence  of  transient  data.  The  flow 
method  usually  does  not  perform  well  in  the  latter  situation  and  often  pro¬ 
vides  a  model  equation  with  unstable  dynamics  instead  of  a  limit  cycle.  The 
advantage  of  the  present  “trajectory  method”,  described  in  detail  in  section 
2,  is  that  it  handles  data  on  the  limit  cycle  very  well  even  when  they  are 
noisy;  the  procedure  generally  produces  a  stable  fit. 

It  is  worth  noting  that  both  the  flow  and  the  trajectory  method  can  also 
be  applied  to  chaotic  time  series  with  generally  good  results.  In  the  present 
paper  we  study  limit  cycle  oscillations,  because  there  are  many  experimental 
systems  in  which  this  kind  of  behavior  appears  [14,15]. 

In  section  3  we  present  some  results,  both  for  simulated  data  and  for 
experimental  data  from  human  limb  movements  [16,12].  Data  from  biological 
systems  are  often  very  noisy,  although  it  still  proves  possible  in  many  cases,  to 
describe  the  trajectories  with  low  dimensional  deterministic  model  equations 
[17].  First  we  study  systems  of  single  oscillators  that  may  be  described  in 


a  two  dimensional  state  space.  Since  the  experimental  data  are  only  on  a 
limit  cycle  and  are  very  noisy,  the  fitted  coefficients  sometimes  depend  quite 
strongly  on  parameters  of  the  fit,  such  as  the  length  of  the  included  dataset. 
Nevertheless  the  trajectory  of  the  given  time  series  is  approximated  quite 
well.  If  additional  information  such  as  transient  data  are  available  we  show 
that  the  reconstructed  coefficients  do  not  change  very  much. 

In  section  3  we  also  study  systems  of  two  coupled  limit  cycle  oscillators, 
which  constitute  important  models  for  biology  in  general  [IS]  and  in  biological 
coordination  in  particular.  There,  self-organization  is  seen  when  different 
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coordination  patterns  arise  as  stable  states  of  the  coupled  nonlinear  dynamics 
[17,19,20].  Up  to  now,  however,  the  reconstruction  is  shown  to  be  successful 
for  simulated  data  with  moderate  additional  noise. 

Section  4  is  a  short  conclusion  that  raises  future  issues. 

2.  Reconstruction  from  time  series  and  the  trajectory  method 

As  a  first  step  it  is  necessary  to  chose  a  state  space  for  the  time  series  to  be 
described.  The  dimension  of  this  state  space  has  to  be  at  least  as  large  as 
the  dimension  of  the  given  time  series.  Various  methods  exist  to  calculate 
the  dimension  of  an  experimental  time  series  [21]. 

In  the  systems  investigated  here,  however,  the  dimension  of  the  system 
is  in  most  cases  quite  obvious  (for  actual  calculation  see  Kay,  Saltzman  & 
Kelso  [12])  and  the  physical  meaning  of  the  measured  time  series  gives  rise 
to  the  dimension  of  the  used  state  space. 

The  state  space  can  be  constructed  either  with  a  time  delay  [13]  or  with 
derivatives.  Differentiation  increases  noise,  but  this  state  space  was  used 
to  describe  the  investigated  low  dimensional  systems  because  only  the  first 
derivative  in  the  trajectory  method  and  the  second  derivative  in  the  flow 
method  are  necessary.  The  advantage  of  using  derivatives  to  define  the  state 
space  is  that  the  components  have  physical  meaning.  In  many  cases  also 
the  derivatives  can  be  measured  directly,  thus  reducing  noise.  Moreover,  the 
state  space  components  must  be  uncorrelated,  which  is  often  a  problem  in 
delay-constructed  state  spaces  [8]  if  the  delay  is  not  chosen  appropriately. 

Many  model  equations  of  nonlinear  oscillatory  systems  use  a  polynomial 
series  as  a  favorable  ansatz  [7].  The  flow  vector  field  is  represented  by  polyno¬ 
mials  of  the  state  vector  components  and  the  coefficients  of  these  polynomials 
are  fitted.  Thus,  given  the  d-dimensional  state  space  vector 


X  = 


/ 

Xi 

\xW 


(2.1) 


the  z-th  component  of  the  flow  is  given  by 

X.  =  £  c.,...„xV  (2.2) 

ll  ,l2...-.>d=0 

where  m  is  the  highest  order  of  the  polynomials  and  are  the  fitted 

coefficients. 

Although  the  method  described  in  [7]  fits  the  flow  directly,  in  the  tra¬ 
jectory  method  the  information  included  in  the  trajectory,  the  evolution  of 
the  state  vector  in  time,  is  used.  Starting  from  appropriately  chosen  initial 
states  Xe  {tj),  the  model  equation  is  u.sed  to  obtain  an  estimation  of  the  states 
for  later  times  Xm  {tj  -I-  A</).  The  coefficients  are  fitted  by  minimizing  the 
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distances  between  the  states  predicted  by  the  model  and  the  experimental 
states.  The  quality  function  Q  is  defined  as 

jmas  Jmax 

Q=  L  +  + Ai,)||  (2.3) 

j./=i 

where  jj  ||  is  the  Euclidian  norm,  Xe  are  the  experimental  states,  tj  are 
the  jmax  times  when  the  initial  states  are  taken  and  At/  are  the  Imax  times 
cifter  which  experimental  and  model  states  are  compared,  is  a  function 
of  the  coefficients  c/,. which  are  determined  by  the  minimum  of  Q  as  a 
function  of  the  coefficients: 


Qmin  =  min  Q  (2.4) 

The  difference  between  this  new  trajectory  method  and  the  recent  works 
of  Cremers  and  Hubler  [7]  and  Crutchfield  and  McNamara  [8]  is  that  it 
employs  both  short  time  and  long  time  behavior  for  the  fit.  Experimental 
and  model  trajectory  are  compared  at  Imax  different  times  Af/.  The  “flow 
method”  uses  only  one  small  time  step  or  only  short  time  behavior,  while 
the  method  by  Crutchfield  and  McNamara  uses  only  the  long  time  behavior 
by  using  one  large  time  step. 

In  eq.(2.3)  the  differences  for  all  state  space  components  are  included.  In 
most  cases,  however,  only  the  squares  of  the  differences  between  the  measured 
experimental  data  and  the  corresponding  model  data 

jmas  Jmas 

Q=  (x,„(f,  +  At,)-x,(fi  + At,))'  (2.5) 

j,i=i 

were  used  to  define  the  quality  function.  In  the  experimental  systems 
studied  here,  usually  the  first  derivatives  of  the  data  are  the  other  components 
of  the  state  space,  but  these  were  not  included  in  the  minimized  quality 
function.  For  such  experimental  data,  the  dynamics  of  the  modelequation 
occasionally  became  unstable  if  the  quality  function  Q  (eq.2.3)  was  used, 
because  the  velocities  were  too  noisy  due  to  the  differentiation  process.  On 
the  other  hand,  when  simulated  data  with  a  small  amount  of  noise  were  used 
to  test  the  present  trajectory  method,  distances  in  the  derivatives  according 
to  eq.(2.3)  were  also  taken  into  account,  though  the  quality  of  the  fit  was 
not  changed.  The  goal  of  the  method  is,  of  course,  to  reconstruct  only  the 
experimental  time  series  itself.  It  is  possible  to  include  more  than  one  time 
series  in  the  fitting  process,  e.g.  in  the  coupled  oscillator  case.  Also,  if  the 
experimental  set  up  allows  the  derivative  to  be  obtained  directly,  it  might  be 
favorable  to  also  use  these  data. 

In  order  to  solve  the  minimization  problem  we  used  a  routine  provided  by 
IMSL  [22].  The  calculation  needs  an  initial  model  equation,  i.e.  initial  coef¬ 
ficients  of  the  model  equation,  which  in  most  cases  are  chosen  to  be  without 
any  force.  If  the  maximal  Af/^^^  is  chosen  very  large,  say  about  one  oscilla¬ 
tion  period,  another  first  guess  for  the  initial  coefficients  is  often  necessary. 
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Otherwise  problems  with  divergences  appear.  As  a  first  guess,  model  equa¬ 
tions  that  were  obtained  using  smaller  were  used.  A  disadvantage 

of  this  method  is  the  large  amount  of  computing  time  necessary  for  large 
At;.  Therefore  only  in  cases  in  which  the  flow  method  fails  is  this  method 
favorable. 

The  At;  were  chosen  according  to 

At,  =  t2‘-^  (2.6) 

where  r  is  the  sampling  rate  and  l^ax  was  chosen  in  the  range  from  4  (8 
samples)  to  9  (256  samples).  But  also  other  definitions  of  At,  are  possible, 
for  example  we  used 


At;  =  2/r  (2.7) 

Nevertheless,  no  significant  differences  in  the  quality  of  the  results  were 
found. 

The  initial  states  Xe(tj)  should  cover  the  whole  limit  cycle.  However,  not 
all  experimental  data  can  be  used  as  initial  states  because  the  calculation 
time  would  be  too  large.  Thus  a  delay  of  b  samples  between  two  initial  states 
was  employed. 

Since  the  trajectory  method  was  applied  only  to  single  oscillators  and  to 
systems  of  two  coupled  oscillators,  we  restrict  the  problem  to  these  systems. 
In  the  case  of  the  single  oscillator,  the  problem  reduces  to 

m 

X  =  ^  (2.8) 

»1  .*2  =0 

and  for  m  =  3  only  10  coefficients  need  be  fitted. 

Instead  of  integrating  a  model  differential  equation,  it  is  also  possible  to 
use  a  map,  which  is  much  less  time  consuming.  If  very  small  time  steps  r, 
smaller  than  1%  of  the  oscillation  period  are  chosen,  the  model  equations  for 
a  single  oscillator  are  given  by 


^n+l  —  Xn  ~T  (Xn  d"  X^+l)  (2-9) 

m 

Xn+l  =  X„  OijijX^Xji^ 

•1  .»2=0 

The  approximation  for  x„+i  is  only  valid  for  small  time  steps,  however. 
For  larger  time  steps  the  fit  of  a  model  equation  for  Xn+i  is  also  necessary, 
thereby  eliminating  the  advantage  of  reduced  calculation  time.  Furthermore 
the  model  equations  become  more  complicated,  because  more  coefficients  are 
necessary. 

The  coefficients  for  an  ODE  can  simply  be  estimated  from  the  coefficients 
for  this  map  according  to 
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c«-  (2.10) 

T 

if  the  time  step  is  small  enough. 

In  the  case  of  two  coupled  oscillators  the  model  equations  are  given  by 


X2 


m 

^2,«i »2 «3 «4 i-i  ^2  i2 

<1  .<2.*3.»4=0 


and  for  a  fit  up  to  third  order  70  coefficients  may  be  fitted. 


(2.11) 


3.  Examples  for  the  reconstruction  of  dynamical  systems 

In  this  section  we  demonstrate  the  reconstruction  methods  on  simulated 
model  equations  as  well  as  on  real  experimental  data.  Some  of  the  latter 
were  previously  reported  rhythmic  movements  of  the  single  hand  about  the 
wrist  [12]  and  a  very  lengthy  single  hand  “control”  time  series  obtained  in 
a  study  of  mode-locking  by  Kelso  and  DeGuzman  [23].  Rhythmic  behav¬ 
ior  seems  a  logical  place  to  start  with  the  present  approach  because  it  is 
biologically  meaningful  in  many  species  [14,20]  . 

3.1  Single  oscillator,  simulated  time  series 

Simulated  data  for  the  well  known  van  der  Pol  oscillator  with  the  ODE 

Xg  =  if  —  XgXg  —  Xg  (3.1) 

were  studied  as  an  example  for  a  single  oscillator  with  a  limit  cycle.  The 
trajectory  method  was  applied  both  to  time  series  with  and  without  transient 
data.  Furthermore,  an  additional  gaussian  noise  force  weis  included  in  eq. 
(3.1)  and  again  data  with  and  without  transient  were  fitted.  In  fig.(l  a)  the 
used  trajectory  with  noise  is  shown.  In  the  following  examples  the  order  of 
the  fitted  polynomials  was  always  3. 

For  transient  data  without  noise  the  coefficients  were  reproduced  very 
well,  with  differences  less  than  1%  relative  to  the  coefficients  of  the  simulated 
differential  equation.  The  coefficients  were  not  reproduced  quite  as  exactly 
for  data  that  lacked  transients,  but  the  differences  were  still  quite  small.  The 
results  of  the  fit  both  for  model  ODE  and  model  maps  are  listed  in  Table  1. 

The  reason  for  the  underestimation  of  the  linear  term  in  z  and  the  over¬ 
estimation  of  the  linear  term  in  z  is  not  yet  understood.  One  possibility  is 
that  the  trajectory  method  may  systematically  overestimate  the  stability  of 
a  limit  cycle,  if  the  data  do  not  include  enough  information  about  transient 
behavior.  On  the  other  hand,  in  many  cases  of  fitted  experimental  and  sim¬ 
ulated  data  that  lie  on  a  limit  cycle,  the  flow  method  produces  a  model  with 
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unstable  dynamics.  This  was  the  reason  for  the  development  of  the  trajectory 
method  in  the  first  place.  Generally  speaking,  only  very  little  information 
about  the  stabilizing  terms  is  present  on  the  limit  cycle  itself.  Therefore,  in 
the  absence  of  transient  data  the  method  used  often  influences  the  stability 
of  the  fitted  model. 

The  fit  of  the  data  with  added  noise  also  reproduces  the  coefficients  quite 
well,  if  transient  data  are  included,  see  Table  2.  For  a  time  series  without 
transient  data  the  coefficients  of  eq.(3.1)  were  not  found  and  for  different 
fit  parameters  the  coefficients  were  also  very  different,  see  Table  3.  But  the 
trajectory  of  these  model  equations  was  always  in  good  agreement  with  the 
trajectory  of  the  van  der  Pol  oscillator.  Fig.l  (b)  displays  the  trajectory  for 
the  model  equation  of  Table  2  with  Imax  =  7.  These  results  suggest  that  in 
many  cases  different  model  equations  may  be  valid  for  the  reconstruction  of  a 
given  trajectory.  However,  if  the  information  in  the  time  series  is  larger,  e.g., 
if  transient  information  is  provided,  the  variety  of  possible  model  equations 
can  be  drastically  reduced  and  model  coefficients  obtained  with  different  fit 
parameters  are  more  or  less  equal. 

3,2  Single  oscillator,  experimental  time  series 

The  foregoing  conclusion  applies  also  to  experimental  data  obtained  in  stud¬ 
ies  of  human  rhythmical  movement.  As  will  be  shown,  the  model  coefficients 
cam  also  be  a  function  of  the  length  of  the  time  series  included  in  the  fit. 

One  time  series  was  a  3Hz  motion  collected  at  a  sampling  frequency  of 
200Hz.  The  movement  is  very  smooth  and  the  limit  cycle  is  not  too  noisy  (see 
fig.  2  a).  In  this  case  the  fit  provided  two  different  model  equations,  whose 
coefficients  nearly  always  matched  one  of  these  equations.  Also  the  standard 
deviations  of  the  coefficients  calculated  with  different  fit  parameters  were 
fairly  small.  The  first  equation  Wcis  found  if  only  five  oscillation  periods  were 
included  in  the  data,  while  the  other  equation  was  found  if  more  or  less  the 
whole  time  series  of  about  ten  periods  was  used.  The  two  model  equations 
are  listed  in  Table  4;  simulated  trajectories  are  shown  in  fig. 2  (b)  and  (c). 
Only  a  fit  with  Ati  chosen  according  to  eq.(2.7)  gave  a  model  equation,  which 
was  again  different  to  the  two  model  equations  listed  in  Table  4.  All  these 
equations  describe  the  given  trajectory  very  well  and  not  one  of  them  can 
be  considered  to  be  better  than  the  others.  The  rather  high  noise  level  often 
present  in  the  experimental  data  make  a  decision  about  which  model  equation 
is  favorable  very  difficult  or  even  impossible. 

Also  for  these  experimental  data  the  coefficients  changed  slightly  and 
homogeneously  with  the  maximal  length  of  Afj.  As  mentioned  earlier,  this 
effect  might  be  due  to  the  reconstruction  method  itself.  More  theoretical 
work  is  required  on  this  issue.  The  flow  method  also  gave  a  model  equation 
(see  Table  4)  with  a  stable  dynamics.  However,  the  time  for  the  relaxation 
to  the  limit  cycle  is  larger  than  for  the  trajectory  method  (see  fig.  2  d).  The 
reason  for  the  difference  in  relaxation  time  lies  in  the  different  fit  methods. 

For  experimental  obtained  at  lower  frequencies  (e.g.  1..5  Hz)  the  fit  be- 
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comes  more  problematic.  In  this  case  the  subject  pauses  at  maximal  displace¬ 
ment  (see  fig.3a).  Such  pauses,  which  are  dependent  upon  the  frequency  of 
motion,  cannot  be  modeled  with  a  two  dimensional  deterministic  equation. 
As  mentioned  in  the  introduction  the  mapping  of  the  time  series  into  state 
space  must  be  bijective,  which  is  not  the  case  for  this  time  series.  If  one 
wants  to  also  fit  this  trajectory  another  approach  is  necessary.  Instead  of 
fitting  one  ODE  to  this  trajectory,  it  is  possible  to  fit  the  movements  in 
each  direction  independently  with  one  ODE  for  each  direction.  The  whole 
dynamics  is  then  described  with  two  discrete  movements  and  a  “trigger”, 
which  switches  according  to  the  performed  frequency  from  one  ODE  to  the 
other.  Of  course,  we  are  not  inferring  from  this  description  alone  that  this  is 
the  way  the  nervous  system  performs  the  task. 

The  fit  of  these  discrete  movements  is  not  possible  with  the  trajectory 
method  as  described  in  section  2,  but  a  fit  with  the  flow  method  is  possible. 
However,  the  model  equations  are  only  valid  in  the  corresponding  part  of 
the  state  space.  If  the  initial  conditions  are  not  chosen  appropriately,  the 
trajectory  may  not  relax  to  the  limit  cycle.  Furthermore,  the  force  must  be 
switched  off  when  the  velocity  passes  zero,  and  the  other  ODE  switched  on 
sometime  later.  In  fig.(3b,c)  a  trajectory  of  a  model  equation  is  shown,  both 
as  a  time  series  and  in  state  space.  The  trajectory  is  reconstructed  quite 
well.  A  fit  for  different  parts  of  the  very  long  experimental  time  series  was 
calculated.  However,  the  coefficients  differ,  depending  on  the  part  of  the  time 
series  that  was  used  for  reconstruction. 

3.3  Reconstruction  of  coupled  oscillators 

For  systems  of  two  coupled  limit  cycle  oscillators  only  simulated  data  have 
been  reconstructed  up  to  now.  To  understand  the  difficulties  in  the  recon¬ 
struction  of  this  system  we  first  review  briefly  some  results  from  studies  of 
biological  movements. 

In  the  original  experiments  subjects  were  asked  to  perform  oscillatory 
motions  with  the  two  index  fingers.  Only  two  phase-locked  patterns  prove  to 
be  stable:  one  in  which  the  fingers  move  symmetrically  (homologous  muscles 
contracting  in-phase)  and  the  other  in  which  the  fingers  alternate  (homolo¬ 
gous  muscles  contracting  anti-phase).  The  discovery,  later  examined  in  great 
detail  and  shown  to  take  the  form  of  a  nonequilibrium  phase  transition  (e.g. 
enhanced  fluctuations,  critical  slowing  down  of  the  order  parameter  were  all 
observed)  was  that  the  anti-phase  pattern  lost  stability  at  a  critical  move¬ 
ment  frequency  and  shifted  to  the  more  stable  in-phcise,  symmetric  pattern 
(see  [20]  for  reviews).  Using  concepts  of  synergetics,  the  order  parameter 
dynamics  for  relative  phase  were  identified  and  later  derived  by  nonlinear 
coupling  [17]  The  significance  of  coupled  nonlinear  oscillators  for  biology  in 
general  [18]  and  the  foregoing  results  showing  that  the  stability  and  change  of 
coordinated  movement  behavior  may  be  understood  as  nonequilibrium  phase 
transitions,  make  systems  of  nonlinearly  coupled  oscillators  a  natural  choice 
of  study  for  the  present  reconstruction  methods.  We  limit  our  analysis  to 
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simulated  data  for  the  moment. 

For  a  system  of  equal  oscillators  in  a  phase  locked  state  the  positions 
and  the  velocities  of  the  oscillators  are  strictly  correlated,  even  equal  for  the 
in-phase  movement.  Since  it  is  necessary  to  have  uncorrelated  state  vector 
components,  a  fit  with  the  state  vector  x  =  (xi,  ij,  X2,  ^2)  is  not  possible 
with  the  present  reconstruction  methods.  However,  since  the  positions  and 
velocities  are  strictly  correlated,  a  model  equation  assuming  two  independent 
oscillators  describes  the  given  time  series  as  well  as  a  system  of  coupled 
oscillators.  Such  a  model  may  even  be  considered  more  favorable,  since  many 
fewer  coefficients  are  necessary  to  model  single  oscillators  without  coupling 
than  a  full  system  of  coupled  oscillators. 

To  reconstruct  the  coupling  terms  transient  data  must  be  available.  If  the 
system  is  attracted  to  a  phase-locked  state  that  is  in-phase,  the  state  vector 

^  Xi 

X  = 

Xi  —  X2 
Xi  —  X2 

is  appropriate,  since  the  third  and  fourth  components  are  attracted  to 
zero  and  the  transient  to  the  phase-locked  state  is  easily  visualized.  The 
differential  equation  of  one  studied  example  is  given  by 


(3.2) 


xi  =  — 2. Oxi -f  0.2xi  —  0.4xiXi -t- 0.1x2 -f  0.3x2  (3.3) 

X2  =  — 2.0x2 -f  0.2x2  —  O.4X2X2 -I- O.lxi -I- 0.3xi 

(3.4) 

An  example  of  a  simulated  trajectory  with  additional  noise  is  shown  in 
fig. 4a  and  b  together  with  the  trajectory  of  the  fitted  model  with  the  same 
initial  conditions.  These  simulations  may  be  compared  with  the  simulations 
in  [17]  and  especially  [24],  which  do  not,  of  course,  reconstruct  the  full  tra¬ 
jectory.  For  the  fit,  both  the  flow  and  the  trajectory  method  were  used  and 
both  reconstruct  the  trajectory  very  well.  Notably,  in  systems  with  a  shorter 
transient  a  good  reconstruction  is  also  possible.  Whether  a  reconstruction 
is  possible  or  not  depends  again  on  the  length  of  the  available  transient  and 
the  noise  in  the  experimental  data. 

The  coefficients  of  the  simulated  differential  equation  are  only  recon¬ 
structed  when  the  additional  noise  is  smaller  than  in  the  plotted  example. 
A  related  problem  appears  when  the  model  is  simulated  using  other  initial 
conditions.  If  the  simulated  time  series  is  too  noisy  or  the  transient  too  short, 
the  model  equation  will  often  produce  an  unstable  trajectory.  The  reason 
is  that  the  flow  vector  field  is  reconstructed  well  in  the  subspace  for  which 
experimental  or  simulated  data  are  available,  but  poorly  outside  this  region. 
For  the  systems  of  coupled  limit  cycle  oscillators  studied  here,  only  data  on 
the  limit  cycle  are  available,  and  the  model  equation  is  valid  only  close  to  the 
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limit  cycle.  In  two  dimensional  systems  with  a  single  oscillator  this  problem 
is  not  so  important,  although  divergences  appear  when  the  initial  conditions 
are  set  further  away  fom  the  limit  cycle.  In  higher  dimensional  systems  such 
issues  are  much  more  important. 

Finally  we  want  to  mention  that  a  reconstruction  of  the  dynamics  of 
coupled  oscillators  is  possible  also  for  other  systems  that  are  not  phase- 
locked.  There  the  results  are  often  better  since  a  larger  subspace  of  the  state 
space  is  reached. 

4.  Conclusion 

In  a  number  of  cases  (both  simulated  and  experimental)  we  have  shown  that 
it  is  possible  to  reconstruct  the  dynamics  of  a  system  from  an  experimental 
time  series  using  a  novel  method,  the  trajectory  method.  With  reconstructing 
methods  the  information  included  in  the  time  series  is  reduced  to  a  very  small 
number  of  coefficients  of  an  ODE,  including  the  same  information  as  the 
time  series.  On  the  other  hand,  we  have  drawn  attention  to  the  drawbacks 
of  such  a  reconstruction  technique  2is  well  as  its  advantages  in  comparison  to 
other  methods  in  literature.  A  number  of  questions  remain  that  are  worth 
addressing. 

One  concerns  how  the  amount  of  information  that  is  necessary  to  produce 
a  successful  fit  can  be  estimated  and  measured.  This  is  a  crucial  issue: 
As  mentioned  above,  fitted  model  equations  generally  cannot  provide  the 
relaxation  time  to  the  limit  cycle  if  no  transient  data  are  available.  How 
much  information  must  be  included  in  the  time  series,  so  that  predictions  for 
global  measures  such  cis  the  relaxation  time  are  possible? 

A  further  issue  concerns  the  influence  of  the  used  fit  method  and  the  used 
fit  parameters  on  the  results  both  of  which  need  to  be  studied  systematically. 
Finally,  the  fit  method  should  give  a  minimal  model  equation  for  the  experi¬ 
mental  time  series,  where  minimality  is  to  be  measured  in  terms  of  order  and 
dimensionality  of  the  fit  and  the  quality  of  the  description  of  the  time  series. 

5.  Survey 

The  essence  of  complex  biological  systems  that  possess  many  degrees  of  free¬ 
dom  is  that  they  can  temporarily  assemble  their  components  into  a  much 
lower  dimensional  structure  whose  dynamics  are  nonlinear  and  hence  capable 
of  exhibiting  complicated  behavior,  among  which  are  oscillations,  bifurcations, 
multistability,  even  deterministic  chaos.  Thus  the  brain,  with  neurons 
and  neuronal  connections.  ~  10^  cell  types,  hundreds  of  active  chemicals  is  a 
materially  complex  system,  whose  dynamics  nevertheless  can  often  be  demon¬ 
strated  to  be  much  lower  dimensional  (for  review  see  e.g.  [6]).  Similarly,  the 
motor  system  (ignoring  for  the  moment  its  neural  and  vascular  support  pro¬ 
cesses)  has  minimally  792  muscles  and  100  joints  to  coordinate  during  the 
course  of  typical  activities.  Watch  a  concert  pianist,  for  example,  and  you 
ivill  see  a  continuous  blending  of  posture,  expression,  and  articulated  biman- 
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ual  performanc:’.  Like  music  itself,  the  artist’s  performance  is  coordinated 
on  many  time  scales. 

To  guarantee  stability  of  behavior  on  the  one  hand  and  flexibility  on  the 
other,  the  behavioral  dynamics  must  be  nonlinear.  In  fact,  in  the  last  few 
years  all  the  signature  features  of  nonequilibrium  phase  transitions,  i.e.  where 
patterns  form  in  a  spontaneous,  self-organized  fashion  under  non-specific 
changes  in  control  parameters  [5]  have  been  observed  and  modeled  in  hu¬ 
man  movement  coordination  tasks  (e.g.  [16,17,19,20]).  Once  the  macro¬ 
scopic  behavioral  patterns  and  their  dynamics  have  been  found  (through  the 
methodology  of  phase  transitions  or  qualitative  changes  in  behavior),  they  can 
be  derived  or  synthesized  by  cooperatively  coupling  nonlinear  oscillators  (e.g. 
[17,24]).  Such  cooperative  coupling  affords  both  a  rigorous  definition  of  the 
relation  between  levels  of  description  and  indicates  how  the  compression  of 
degrees  of  freedom  occurs,  e.g.  from  the  limit  cycle  dynamics  of  the  individual 
componer  '  to  the  point  attractor  dynamics  of  the  collective  variable,  or  or¬ 
der  parameter  relative  phase.  Thus,  the  reduction  is  not  to  some  fundamental 
unit  or  level  of  analysis,  but  rather  to  find  laws  at  one  level  and  derive  them 
from  another. 

In  this  regard,  it  is  quite  clear  that  nonlinear  oscillators  play  quite  a  cen¬ 
tral  role  at  many  levels  of  description  and  throughout  the  natural  sciences. 
A  large  variety  of  complex  systems  effectively  reduce  their  degrees  of  free¬ 
dom  to  the  frequency  and  phase  or  mode-locking  dynamics  of  coupled  non¬ 
linear  oscillators  ([14]  for  recent  review).  One  of  us  is  inclined  to  think 
that  this  is  a  strategy  for  the  nervous  system  too  which  must  perform  tasks 
that  require  the  coordination  among  many  spatially  and  temporally  distributed 
processes  (see  e.g.  [23]).  Such  a  coordinative  design  is  not  only  apparent 
in  voluntary  behavior  but  in  many  other  “systems”,  ranging  from  so-called 
locomotory-respiratory  rhythms,  cardiac  rhythms,  cardiovascular-respiratory- 
somatomotor  interactions,  to  individual  neurons  such  as  the  squid  axon  and 
so  forth.  Stationary  states,  oscillations,  bifurcations,  intermitt ency,  multi- 
stabilicy,  chaos  -  in  other  words,  rich  dynamics  that  depend  on  the  region  in 
parameter  space  in  which  the  system  lives  -  are  ubiquitous  features  of  cou¬ 
pled  or  forced  nonlinear  oscillators.  Thus,  in  the  present  article,  the  nonlin¬ 
ear  oscillator,  exemplified  by  data  obtained  from  studies  of  human  rhythmical 
behavior,  is  the  chief  target  of  study. 

Of  course,  this  level-independent  strategy  on  complex  systems  means  that 
one  must  be  able  to  place  the  system  under  study  in  a  context  that  will  en¬ 
able  one  to  deduce  lawful  structure  from  real  data  that  are  often  noisy,  the 
so-called  inverse  problem.  Above,  we  mentioned  using  phase  transitions  or 
behavioral  instabilities  as  a  special  entry  point  where  laws  may  be  found.  In¬ 
stabilities  serve  to  demarcate  behavioral  patterns,  thereby  allowing  a  precise 
identification  of  order  parameters  for  patterns  and  their  nonlinear  dynamics 
[5,19,20].  But  what  other  methods  may  allow  one  to  obtain  dynamical  laws 
(equations  of  motion)  directly  from  measurements?  And  what  steps  can  be 
taken  to  optimize  these  descriptions?  Many  tools  and  techniques  have  been 
employed  in  efforts  to  characterize  the  dynamics  of  complex  .systems,  includ- 


time  series 


13 


ing  dimension,  metric  entropy,  Lyapunov  exponents  and  so  forth  (see  e.g. 
[6]).  Here  we  introduce  a  new  method  that  reconstructs  the  entire  flow  vector 
field  from  the  individual  trajectory  of  the  system,  i.e.  the  evolution  of  a  point 
in  a  low  dimensional  state  space.  The  reconstructed  flow  vector  or  ordinary 
differential  equation  constitutes  a  global  description,  i.e.  is  valid  on  the  whole 
subspace  of  the  state  space  where  the  experimental  data  reside.  We  discuss 
the  advantages  and  disadvantages  of  this  method  in  the  article  itself. 

Reconstruction  of  the  underlying  dynamics,  of  course,  does  not  nece.ssarily 
provide  a  physical  explanation  or  enhance  scientific  insight,  per  se.  Neverthe¬ 
less,  when  complemented  by  theoretically  motivated  experimental  observations 
it  may  still  be  possible  to  connect  model  equations  to  deeper  levels  of  descrip¬ 
tion.  Certainly,  in  the  history  of  science  there  are  examples  where  an  initially 
guessed  mathematical  formalism  and  eventual  physical  interpretation  led  to 
the  same  results.  In  the  case  of  complex  systems,  where  it  is  often  not  possible 
to  derive  a  mathematical  structure  from  first  principles,  these  reconstruction 
methods  could  play  an  important  role. 

One  possible  application  for  methods  of  reconstruction  is  in  the  field  of 
controlling  dynamical  systems.  It  may  be  possible,  for  instance,  to  control 
very  complex  systems  if  the  dynamics  of  the  system  can  be  described  with  a 
small  number  of  order  parameters.  Of  course,  for  such  control  it  is  necessary 
to  describe  the  dynamics  of  the  uncontrolled  system  as  exactly  as  possible. 
Here  reconstruction  can  be  an  important  tool. 

The  goal  of  the  present  work  is  not  simply  to  analyze  experimental  mate¬ 
rial,  but  to  show,  in  general,  what  is  or  might  be  possible  in  the  field  of  recon¬ 
struction  from  time  series.  The  described  methods  can  be  used  in  all  fields  of 
scientific  research,  where  dynamical  variables  are  observed.  Biological  move¬ 
ments,  the  dynamics  of  electrical  potentials  in  nerves,  fluid  mechanics  or  the 
dynamics  of  chemical  reactions  are  some  examples.  The  central  ansatz  of 
the  method  is  simply  that  any  observations  are  produced  by  a  dynamical  sys¬ 
tem  in  the  presence  of  noise.  The  aim  is  to  provide  a  minimal  model  that 
reproduces  observed  behavior.  We  hope  that  the  given  examples  will  help  re¬ 
searchers  to  think  about  possible  applications  of  the  present  methods  in  their 
field  of  interest. 
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modeltype 

^max 

Co,i 

Cl.O 

^2.1 

other  c 

ODE 

4 

1.147 

-0.9910 

-1.041 

<  0.065 

ODE 

5 

1.129 

-0.9729 

-1.028 

<  0.042 

ODE 

6 

1.135 

-0.9564 

-1.024 

<  0.035 

ODE 

7 

1.128 

-0.9472 

-1.021 

<  0.030 

map 

4 

1.199 

-1.102 

-1.157 

<  0.019 

map 

5 

1.015 

-0.9321 

-0.9859 

<  0.030 

map 

6 

1.006 

-0.9296 

-0.9802 

<  0.036 

map 

8 

1.009 

-0.9170 

-0.9781 

<  0.041 

map 

9 

1.013 

-0.9175 

-0.9795 

<  0.039 

Table  1:  The  coefficients  for  the  fit  of  a  simulated  van  der  Pol  oscillator 
(eq.  3.1).  r  was  0.067  for  fit  with  an  ODE  and  0.02  for  map.  The  other 
parameters  were:  6  =  7,  jmax  =  100,  the  At(  were  chosen  according 
to  eq.(2.6).  The  used  data  were  without  noise  and  without  tramsient. 
Map  coefficients  are  rescaled  according  to  eq.(2.10) 


^max 

Cq.I 

^1,0 

^2.1 

other  c 

4 

1.594 

-1.230 

-1.125 

<  0.317 

5 

1.395 

-1.476 

-1.103 

<  0.308 

6 

1.153 

-1.293 

-0.9650 

<  0.194 

7 

0.9140 

-0.8094 

-0.8282 

<  0.045 

Table  2:  The  coefficients  for  the  fit  of  a  simulated  van  der  Pol  oscillator 
(eq.  3.1)  with  an  ODE.  The  time  series  was  with  noise  (see  fig.l  a)  and 
with  transient  data.  The  other  parameters  were:  r  =  0.067,  6  =  7, 
jmax  =  100,  the  Alt  were  chosen  according  to  eq.(2.6). 
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coefficient 

^max  —  4 

imax  —  4 

^max  —  5 

Imax  =  6 

^max  —  ^ 

Oo.o 

0.306 

-0.022 

-0.433 

-0.089 

-0.232 

Oo,l 

4.063 

3.173 

4.349 

4.2.35 

3.814 

Co,2 

-0.117 

-0.034 

0.102 

0.006 

0.043 

Co, 3 

-0.638 

-0.455 

-0.708 

-0.693 

-0.621 

Cl.O 

-1.338 

-1.708 

-0.868 

-0.391 

-0.529 

^1,1 

0.131 

0.009 

-0.058 

-0.045 

0.003 

Cl, 2 

0.990 

0.829 

0.943 

0.750 

0.662 

C2,0 

-0.046 

0.010 

-0.133 

0.061 

0.092 

C2,l 

-2.089 

-1.883 

-2.054 

-1.878 

-1.750 

C3,0 

-0.022 

0.104 

-0.137 

-0.237 

-0.187 

Table  3:  The  coefficients  for  the  fit  of  a  simulated  van  der  Pol  oscillator 
(eq.  3.1),  as  in  table  (2)  but  without  transient  data.  The  second 
column  with  Imax  =  4  is  the  result  of  a  fit  with  jmax  =  200.  Note  the 
large  differences  in  the  linear  coefficient  ci,o. 


coefficient 

short  time  s. 

SD 

long  time  s. 

SD 

2/r 

flow  method 

co,o 

1.227 

0.528 

-1.169 

0.198 

-0.856 

-1.368 

Co,l 

0.916 

0.092 

0.219 

0.729 

0.424 

Co,  2 

-0.599 

0.132 

-0.008 

0.029 

-0.075 

0.055 

Co,3 

-0.189 

0.032 

-0.064 

0.037 

-0.077 

-0.009 

Cl,0 

-4.825 

0.366 

-4.359 

0.246 

-5.524 

-3.856 

Cl,l 

0.849 

0.040 

-0.827 

0.068 

1.019 

0.800 

Cl, 2 

0.063 

0.061 

0.074 

0.065 

0.263 

0.034 

C2,0 

-1.094 

0.399 

0.897 

0.138 

0.875 

0.976 

C2,l 

-1.083 

0.107 

-0.615 

0.1.50 

-0.760 

-0.448 

C3,0 

0.4.52 

0.127 

0.670 

0.147 

1.488 

0.324 

Table  4:  The  coefficients  for  the  fit  of  an  experimental  time  series,  a 
single  hand  movement  as  in  fig.(2)  with  the  trajectory  method.  The 
standard  deviation  was  calculated  with  the  coefficients  from  fits  with 
6  different  fitparameters  in  the  case  of  the  short  time  series  and  9 
different  fitparameters  in  the  case  of  the  long  time  series.  For  these 
fits  the  Ati  were  chosen  according  toeq.(2.6).  Some  of  the  coefficients 
are  small  and  can  be  neglected,  Ci,2  in  both  cases  and  co,2  for  the  long 
time  series.  The  cubic  term  in  the  velocity  co.3  is  small  but  important 
for  the  stability  of  the  model  equations.  Furthermore  the  result  with 
At/  chosen  according  to  eq.(2.7)  is  listed.  For  comparison,  the  result 
of  a  fit  with  the  flow  method  is  listed,  the  difference  is  basically  in  the 
cubic  term  in  the  velocity  co,3,  which  is  very  small. 
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Figure  1:  (a)  A  trajectory  of  a  simulated  van  der  Pol  oscillator  (eq. 
3.1)  with  noise  in  a  state  space  with  position  and  velocity  as  state 
vector  components.  The  trajectory  up  to  the  point  marked  with  t 
was  excluded  for  fits  without  transient  data. 

(b)  A  trajectory  in  the  same  state  space  for  the  model  equation  of 
table  (3)  with  /max  =  7.  Although  the  coefficients  are  very  different 
from  the  coefficients  of  eq.(3.1)  the  shape  of  the  trajectory  simulated 
with  the  model  equation  is  very  similar  to  the  trajectory  of  the  van 
der  Pol  oscillator. 


Figure  2:  (a)  A  trajectory  of  an  experimental  single  hand  movement 
(around  wrist)  at  a  frequency  of  3Hz  in  arbitrary  units  in  a  state  space 
with  position  and  velocity  as  state  vector  components.  Trajectories 
in  the  same  state  space  for  the  model  equations  of  table  (4),  (b)  long 
time  series,  (c)  short  time  series,  (d)  flow  method. 


Figure  3:  (a)  Displacement  for  a  single  finger  movement  at  a  frequency 
of  1.5Hz  in  arbitrary  units.  Note  the  pause  at  maximal  displacement. 

(b)  Time  series  of  the  model  equation  obtained  by  the  flow  method. 

(c)  Same  time  series  as  (b)  in  state  space  with  position  and  velocity 
as  state  vector  components. 


Figure  4;  (a)  (b)  Trajectory  of  eq.{4)  in  the  state  space  given  by 
(2).  The  slightly  distorted  line  is  the  simulated  time  series  with  noise, 
the  straight  line  is  a  simulation  of  the  model  equation  with  the  same 
initial  condition.  For  ease  of  visualization  in  (c)  and  (d)  we  plot  the 
positions  of  both  oscillators  against  each  other  to  show  the  transition 
from  anti-phase  to  in-phase  pattern,  (c)  is  a  simulation  of  eq.(4) 
(without  noise)  and  (d)  is  a  trajectory  of  the  model  equation  obtained 
by  the  trajectory  method. 
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